arXiv:1507.01894v2 [cs.CE] 8Jul2015 


On pore-scale modeling and simulation of reactive transport in 3D 

geometries 


Oleg Zahra Lakdawala®’'^, Katherine H. L. Leonard*^, Yavor Vutov"^ 

‘^Fraunhofer Institute for Industrial Mathematics ITWM, Germany. 

^Numerical Porous Media SRI Center, King Abdullah University of Science and Technology, Kingdom of Saudi Arabia. 

‘^DHI-WASY GmbH, Berlin, Germany. 

‘^Institute of Information and Communication Technologies, Department of Scientific Computations, Bulgarian Academy 

of Sciences, Sofia, Bulgaria. 


Abstract 

Pore-scale modeling and simulation of reactive flow in porous media has a range of diverse applications, 
and poses a number of research challenges. It is known that the morphology of a porous medium has 
significant influence on the local flow rate, which can have a substantial impact on the rate of chemical 
reactions. While there are a large number of papers and software tools dedicated to simulating either 
fluid flow in 3D computerized tomography (CT) images or reactive flow using pore-network models, little 
attention to date has been focused on the pore-scale simulation of sorptive transport in 3D CT images, 
which is the specific focus of this paper. Here we first present an algorithm for the simulation of such 
reactive flows directly on images, which is implemented in a sophisticated software package. We then 
use this software to present numerical results in two resolved geometries, illustrating the importance of 
pore-scale simulation and the flexibility of our software package. 

Keywords: Reactive transport modeling. Pore-scale model. Finite volume method, CFD, Surface 
reactions. Filtration 


1. Introduction 


Understanding and controlling reactive flow in porous media is important for a number of environmen¬ 
tal and industrial applications, including oil recovery, fluid filtration and purification, combustion and 
hydrology. Traditionally, the majority of theoretical and experimental research into transport within 
porous media has been carried out at macroscopic (Darcy) scale. Despite the progress in developing 
devices to perform experimental measurements at the pore-scale, experimental characterization of the 
pore-scale velocity, pressure and solute fields is still a challenging task. Due the influence of the pore- 
scale geometry on the processes of interest, direct numerical simulation (DNS) promises to be a very 
useful computational tool in a wide range of fields, and in combination with experi mental studies, ca n 
be used to determine quantities of interest that are not experimentally quantifiable ([Blunt et al.L 12013 ). 

Significant progress over the past 10-15 years in the pore-scale simulation of single phase flow has 
resulted in the computation of permeability tensors for natural and technical porous media becoming 
a standard procedure. A number of academic as well as commercial software tools, capable of process¬ 
ing 3D CT images in addition to virtually generating porous media, are available, for example Aviso, 
GeoDict and Ingrain. Most of those software tools have the additional ability of simulating two-phase 
immiscible flow at the pore-scale directly on a computational domain obtained through the segmentation 
of 3D CT images, often using the lattice Boltzmann (LB), the level set or volume of fluid methods. In 
contrast, substantially less work on the DNS of reactive flow has been performed, and only a few software 
tools with this capability exist. A limited number of computational studies examinin g reactive transport 


where the reactio ns only occur within the fluid phase (and not at a surface) exist ([Molins et al.L 12012 : 
Shen et l2Qll[ ). In contrast, the literature and computational tools examining full 3D reactive flow 
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where the reactions occur at the pore wall (surface reactions), is sparse. The majority of existing stud¬ 
ies and available numerical simulation packages descr ibing pore-s c ale so r ptive transport are base d on 

and 


les and availaole numerical simulation pacKages descr ioing pore-s c ale so rp 
pore-network mathematical models (see, for example. iRaoof et all ( 2Q1Q[) : H 


IVarloteanx et al. 


literature therein), where the geometry needs to be converted into an idealized series of connected pores 
and throats to represent the porous medium. However, duri ng this proce s s, info rmat ion on the morphol¬ 
ogy o f the underlying media can be lost (see, for example, iRaoof et al. ( 2QlQh and iLichtner and Kang 
(|2QQ7h b In this paper we present an algorithm and a sophisticated software package, called Pore-Chem, 
which uses cell centered finite volume (FV) methods to numerically solve 3D solute transport with sur¬ 
face reactions at the pore-scale. In particular the software package has the ability to solve the systems of 
equations modeling colloidal reactive transport on a geometrical domain obtained directly through imag¬ 
ing techniques, such as computerized tomography, which allows for a very accurate spatial description 
of the computational domain. 

In this paper we describe the transport of a generic solute through the Navier-Stokes (NS) system 
of equations coupled to a convection-diffusion (CD) equation. The CD equation is complemented by 
boundary conditions which describe various types of surface reactions comprising a Robin boundary 
condition for the solute coupled to an ordinary differential equation (ODE) describing the dynamics 
of the adsorption at the pore wall. To represent the computational domain, a voxelized geometry is 
considered, where each individual voxel is either solid or fluid. The solute transport model is chosen 
due to its applicability to a broad range of problems. In addition, our goal is to describe the transport 
and reaction of sub-micron particles, for which inertial effects of the individual particles are negligible. 
Discrete models, where each particle is modeled as an individual entity, are necessary when considering 
larger particles, for example with a radius greater than one micron, for which inertial effects become 
important. Several commercial software packages, for example GeoDict ( GeoDict. 2012-20141) . solve a 
range of discrete mathematical models describing colloid transport and adsorption. However, numerical 
simulation of these models is often significantly more computationally expensive than a continuous 
mathematical model and accounting for different reaction kinetics is usually not possible in such packages. 

Reactive flow in porous media is intrinsically a multiscale problem. The goal of our developments is 
to support problems where scale separation is possible and in cases where it is not possible. The first 
case, where the separation of scales is viable, is usually the focus of asymptotic homogenization theory. 
In the second case, when scale separation is not possible, numerical upscaling methods like multiscale 
finite element methods are often applied. During the homogenization procedure, when applicable, certain 
assumptions are imposed, allowing for the derivation of macroscopic (Darcy scale) equations from the 
microscale formulation, with effective parameters , such as the p ermeability and the effective reaction 
rate, obtained through solution of a cell problem ( HornnnS. 1 19971) . A number of studies have employed 
homogenization theory to derive a macroscopic description of sorptive reactive transport for particular 
parameter regimes. The homogenization of solute transport in porous media in the prese nce of surface 
reactio ns has been performed for both h igh Peclet numbers (convection dominated regim e) (lAllaire et al. 


2QlQblla : Allaire and Hutridurga ^ 20121) . and when the Peclet number is of order one ( Hornnn^ 1994 : 


Kumar et al. ■ 120131: 1 van Duiinl . Il99l[) . In addition to being able to solve cell problems in a number of 


settings, our software has the ability of solving a much broader class of problems at the pore-scale, 
without being restricted by the assumptions required during homogenization. Furthermore, it provides 
the possibility to study various different types of surface reactions described by different kinetics. 

The remainder of the paper is organized as follows. In Section [5] the mathematical model is presented 
and cast into dimensionless variables. The method of achieving a numerical solution to the system of 
equations is outlined in Section [3l and illustrative results using this method are presented in Section 01 
Finally, conclusions are drawn in Section 0 


2. Mathematical model 

We now detail the mathematical model, which describes the transport and reaction of a generic solute 
at a 2D interface within a 3D pore-scale resolved geometry, where we assume that the number of solute 
particles is sufficiently large that representation within a continuum framework is valid. 

Let us denote the spatial domain of interest by ^2, an open subset of We assume that we 
can decompose Q into a solid domain, denoted by and a fluid domain, denoted by such that 
Q = Denoting the external boundary of our domain, being the closure of D, by dQ, we partition 

this into an inlet, an outlet, and external walls, so that 


dQ = dQin U dQout U Qf^wall- 
2 


( 1 ) 













































We note that, although we here consider only one inlet and one outlet, extension to consider multiple 
inlets and outlets is straightforward. Finally, we denote the inter facial boundary between the fluid and 
solid portions of the domain by F = ftff] fig- To allow for different types of reactive boundary conditions 
to be described, we decompose F into N >1 different boundary types as follows 

r = ufL-o^Ti, (2) 

where F^ has distinct properties. Making such a decomposition enables us to attribute different reaction 
rates or kinetic descriptions to different portions of the domain, allowing for different types of solid 
material. 

In order to describe the flow of the water through the membrane, by appealing to the conservation 
of momentum, the incompressible NS equations are used: 

V • V = 0, (3b) 


where v(x, t) and p(x, t) are the velocity and pressure of the fluid respecti vely, while M > 0 and p > 0 
are the viscosity and the density of the fluid which we assume are constants ( Bear! . 1988 : Achesonl . 1995 b 
Suitable boundary conditions on dft are given by 


V = Vin, 

X G BClin, 


(4a) 

P ~ ^out 

X G 

t > 0. 

(4b) 

v = 0, 

X G 


(4c) 


where Vin is the inlet velocity with ||Vin|| > 0, Pout is the pressure at the outlet, and n is the outward 
pointing normal to the boundary dQ. Although here we have used no-slip and no-flux flow conditions 
for the external walls, symmetry or periodic boundary conditions can also be imposed which may be 
more appropriate depending on the problem to be solved. Further boundary conditions are required to 
be specified on the reminder of the boundary to Clf, being the fluid-solid interface. To allow for the slip 
of flow along the fluid-solid interface, and the inclusion of additional effects such as charged fluids or 
matrices, we use the Navier-Maxwell slip conditions, given by 


V • n = 0, V • t = PiU 


(Vv + (Vv)^) ■ t, 


X G Fi, t > 0, 


( 5 ) 


where /3i is the slip length onx G F^ fori = 0,1,...A^ — 1 measured per unit leng th, t is any unit tan gent 
to the surface such that t • n = 0, and the superscript T denotes the transpose (|Langa et ah . 2007). In 
the case that A = 0 for some i, then the standard no-slip and no-flux boundary conditions for the flow 
are enforced on F^. We specify initial conditions through 


v(x,0) = vo(x), p(x,0) =po(x), xGf^/, (6) 

where vq and po are known functions. We discuss the choice for these in Section [3l 

We denote the concentration of the solute within the fluid by c(x, t), measured in particle number 
per unit volume. Appealing to the conservation of mass and assuming no fluid-phase reactions occur, 
the spatio-temporal evolution of the solute concentration is given by 

Be 

— + V • (vc) = DV • (Vc), xGf^/,t>0, (7) 

where D > 0 is the solute diffusion coefficient which we assume to be scalar and constant. We assume a 
known concentration of the solute at the inlet, and prescribe zero flux of the solute at the outlet and on 
the external walls as follows: 


C — Cin, 

Vc • n = 0, 


X G BClin, (8a) 

t > 0, ^ ^ 

X G BClout u (8b) 


where cin > 0 is assumed to be constant. 
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2.1. Models for surface reactions 

We are required to specify boundary conditions for c(x, t) on x G to describe the surface reactions 
occurring here. In general, there are two stages of adsorption of a particle from the bulk solution to the 
solid surface. The first stage involves the diffusion of particles from the bulk solution to the subsurface 
and the second stage then involves the transfer of particles from the subsurface to the surface. After the 
adsorption of a molecule at the interface, there i s a reorientation of the c olloid molecules at the surface, 
which results in a change in the surface tension ( Kralchevskv et aP. I2QQ8I) . Assuming that both the rate 
of diffusion of the particle from the bulk to the subsurface, and the rate of the transfer of the particles 
from the subsurface to the surface are important in determining the overall rate of reaction, we use a 
mixed kinetic-diffusion adsorption description, given by 


D\/c • n = 


dm 

dt 




X e Tj. t > 0, 


( 9 ) 


Here is the surface concentration of the particle under consideration (iKralchevskv et aill2flfl8 h 


measured in units of number per unit surface area, which contrasts with c(x, t) being measured in units 
of number per unit volume. The function fi{c^m) describes the kin etics of the rate of change of the 
surface concentration on the reactive boundary for i = 0,... A" — 1 ( Panov et al.L [20021) . Equation ([9]) 
states that the change in the surface concentration is equal to the flux across the surface, where the 
movement from the bulk to the surface is termed adsorption^ and movement from the surface to the bulk 
is termed desorption. If T^ is nonreactive, for some i, then /^ = 0, so the adsorbed concentration on 
this boundary type remains constant, and a no-flux boundary condition for the solute concentration is 
prescribed. For reactive boundaries, the choice of fi and its dependence on c and m is highly influential 
in correctly describing the reaction dynamics at the solid-fluid interface. A number of different isotherms 
exist for describing these dynamics, dependent on the solute attributes, the order of the reaction, and 
the interface type. 

The simplest of these is the Henry isotherm, which assumes a linear relationship between the surface 
pressure and the number of adsorbed particles, and takes the form 


fi = i^a,iC - nd,im, 


X G Fi, t > 0. 


( 10 ) 


Here Ka,i > 0 is the rate of adsorption, measured in unit length per unit time, and Kd,i > 0 is the rate of 
desorption, measured per unit time, at reactive boundary type i for i = 0,... A — 1. Equation (uni states 
that the rate of adsorption is proportional to the concentration of particles in solution at the reactive 
surface. As such, the rate of adsorption predicted does not saturate at higher surface concentrations. 
However, physically we expect the rate of adsorption to decrease as the quantity of adsorbed particles 
increases and the available surface area for adsorption decreases. Even though the Henry isotherm 
predicts no limit to surface concentration and does not model any interaction between the particles, it 
has been used in a large number of analytical studies due to its linearity. 

The Langmuir adsorption isotherm was the first to be derived mathematically, and is suitable to 
describe the adsorption of a monolayer of localized non-ionic non-i nteracting m olecules at a 2D solid 
interface, and a derivation from statistical physics may be found in ( Baretl [IqGQI) . It is also frequently 
used to describe the adsorption of molecules at a solid-liquid interface, and is described by 


fi — 



i^d,im, 


X G Fi, t > 0. 


( 11 ) 


Here moo,i > 0 is the maximal possible adsorbed surface concentration, measured in number per unit 
area, at reactive boundary type i for i = 0,... A — 1. In comparison to the Henry isotherm, the 
Langmuir isotherm predicts a decrease in the rate of adsorption as the adsorbed concentration increases 
due to the reduction in available adsorption surface. The Henry isotherm, given in Equation , is a 
linearization of Equation (HU, explaining why it produces an accurate representation only at low surface 
concentrations. 

More complex descriptions of the reaction kinetics exist to describe non-localized adsorption and 
particles which interact. For example, the Frumkin isotherm describes localized adsorption where particle 
interaction is accounted through an additional parameter: 


fi — l^a,i^ 



Kd,im exp 
4 



X G Fi, t > 0. 


( 12 ) 





























where k is the Boltzmann constant and T is the temperature in Kelvin. The parameter /3 > 0 describes 
how cooperative the reaction is, and is related to the interaction energy between the particles. In the 
case that /d = 0, the Langmuir isotherm is recovered. The Frumkin isotherm is infrequently used in the 
mathematical modeling of colloidal transport, which may be due its nonlinearity and the difficulties in 
determining the interaction energy between the particles. 

We make the assumption that the adsorption or desorption of our solute does not alter the position 
of the reactive boundary, which in the case of small volumes of particles being adsorbed is sufficiently 
accurate. By the conservation of mass, such an assumption implies any adsorption or desorption on the 
surface is represented by a corresponding increase or decrease in the density of the solid material through 
time. In some cases, for example when the molecules are big or the number being adsorbed is large, 
i nterface evolution needs to be considere d an d may be done in a simila r manner to iTartakovskv et ahl 
(2008); Roubinet and Tartakovskv ( 20131) and Boso and Battiat^ ( 20131 ). This is particularly important 
in the application of rock dissolution and precipitation, where large geometrical changes are observed. 

To close the system of equations, we impose the initial conditions 


c(x, 0) = co(x), X G O/, m(x, 0) = mo,i(x), x G T^, (13) 

where cq and mo,i are known for i = 0,1,... A/" — 1. 

Our problem is, therefore, described by two systems of equations with one-way coupling; the incom¬ 
pressible NS equations, described by ® - ® and the CD equation described by © - with reactive 
boundaries conditions m, initial conditions m, and a description of the reaction kinetics, for example 
Equation (nni), (CH) or m- We now cast the equations into dimensionless variables, before detailing the 
methods used to obtain a numerical solution. 


2.2. Nondimensionalization 

Using a caret notation to distinguish the dimensionless variable from its dimensional equivalent, we 
let 

X = Lx, V = 'l^nV, P = Pont + pV^p, C = CinC, m = CinLm, M = CinL^M, 

where L > 0 is a typical length scale of the computational domain and Un = ll^inll- As our computational 
domain consists of voxels, the relationship between each voxel and its material property is conserved 
upon nondimensionalization, while the length, area and volume of each voxel scales with L, and 
respectively. Given this, we let U, (ig and U/, with boundaries 9Uin, and T^ for i = 

0,1,... A/" — 1 represent the dimensionless versions of the equivalent dimensional domains and boundaries, 
where the voxel size is scaled accordingly. In dimensionless variables, we, therefore, have 


dv 


dt 

dc 


1 


— + V • Vv = -Vp + VW 


Re 


V-v = 0. 


X G O/, t > 0, 


^, + V^(vc) = -V^(Vc), 


(14a) 

(14b) 

(14c) 


where 


Re = 


LpVin 



(15) 


are the global Reynolds and Peclet numbers respectively, being the ratio between the inertial and viscous 
forces and the ratio between advective and diffusive transport rates respectively. Boundary conditions 
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are given by 


V 

= n. 

X G 

a^in, 

(16a) 

p 

= 0 and Vv • n = 0 

X G 


(16b) 

V 

= 0, 

X G 

^^wall 1 

(16c) 

V • n 

= 0, and V • t = /3in • ^Vv + ^ 

X G 

fi, 

i > 0. (16d) 

1 ^ 

drh - , 




-vs n 

= -^ = fiic,m), 

at 

X G 

fi, 

(16e) 

c 

= 1, 

X G 


(16f) 

Vc • n 

= 0, 

X G 

^^out U 

(16g) 


with Pi = In the case of the Henry isotherm, we have 


fi = Da^^^c - Da^^^m, x G t > 0, 
whereas, nondimensionalization of the Langmuir and Frumkin isotherms yields 

m 


fi = 1 


Da^^^m, X G Fi, t > 0, 


and 


fi = DaVc 1 - 


m 


ruo 


— Da^ exp 


(- 


X e Tj, t > 0, 


(17) 


(18) 


(19) 


respectively, for i = 0,1,... A/" — 1, where P = and moo,i = 


numbers are given by 


Da^ 7' — 




kT 

and Da^^ i = 


C’lYiL 

1^’ 


. In ([171) - (dH) the Damkohler 

( 20 ) 


and describe, for each i, the ratio of the rate of reaction (either adsorptive or desorptive) to the rate of 
advective transport. The initial conditions are given by 


v(x, 0) = Vo (x), p{x, 0) = Po (x), c(x, 0) = Co (x), 

m(x,0) = mo,i(x), 


X e 0/, 
X e f j, 


(21a) 

(21b) 


, Vo(x) Po{x)-Pout .... Co(x) TOo,i(x) • ni AT 1 

where vo(x) = —-—, po{x) = -—^-, co(x) =-and mo,i(x) = -— for ^ = 0,1,... -1. 

Idn P^in Qn-L 

We now proceed to discuss the numerical methods used to obtain an approximate solution to our 
dimensionless system of equations given by m, with boundary conditions given by (HH) - (HU) and 
initial conditions specified through m- 


3. Numerical methods 

The full system of equations cannot be solved using analytical techniques, and so numerical methods 
need to be employed to calculate an approximate solution. We employ FV methods to numerically 
solve our system of equations, motivated by their local mass conserving properties. Other methods, for 
example LB or finite difference methods may also be used for solving the flow problem. We note that 
our CD solver is completely compatible with LB methods (the compatibility of our solver with finite 
difference methods depends on the grid selection). Although the authors are not aware of a detailed 
comparison of the performance of FV and LB methods in solving the NS equations at the pore-scale, 
some incomplete internal studies indicate that LB methods can be advantageous for geometries with a 
very low porosity and a high tortuosity, while FV methods are favorable in other cases. 

Due to the one-way coupling between our two systems of equation, the velocity and pressure solutions 
are at steady state. For the sake of generality, we consider the unsteady equations, and begin by solving 
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the system of equations describing the fluid flow, namely ()14ap . (jl4bp along with ()16ap - (|16dp . to obtain 

dir 

a steady state numerical solution where —r = 0 is satisfied. This is achieved using a Chorin fractional 

timestepping method and we refer the reaSr tn icie.is et ai.l ^ and iT.akdawaial ^ for full details 
and for further references, where the methods used are described. 

Once the solution of v is obtained, we proceed to solve the system of equations describing the solute 
transport and reaction, (jl4cp along with the boundary conditions dni - (HI and initial conditions 
m, using a FV method with a cell-centered grid. For the sake of brevi ty, full details of the numerical 
method employed are not given, but we refer the reader to, for example ICauson et al.l ( 2011 ). for more 
detailed information on FV methods. Firstly, dimensionless time is uniformly partitioned into Q time 
points, denoted by with = k (At), where At is the dimensionless timestep size. Then 

the spatial domain, 11, is split into P non-overlapping cubic finite volumes, S/ for / = 0,1,... Pi, which 
span the 3D computational domain, such that Cl = Considering a single representative finite 

volume, S/, we denote its six faces by Pij with center Xj where the subscript j = s^t^b denotes 

the east, west, north, south, top and bottom faces respectively. Integration of (jl4cp over the control 
volume, S/, and time interval upon application of the divergence theorem, yields 


£fc + l 

[ c(x,?+^)dV- [ c(x,?)dV-|- [ [ ■ 

Jbi Jbi JP JdBi 


VC • n dS dr = 


fk+i 

f [ 

Pe Jik Jq]s^ 


Vc • n dV dr. 


where dBi is the boundary of S/, so that dBi = w n st • Denoting the center of the finite 

volume by Xc and using the approximations / (/)(x, t) dP = Aj(j){i^j , i) and / 0(x, t) dP = |S/ |0(xc, t), 

^ Jbi 

for some scalar function cj) for j = n, 5, e, re, t, 6, where Aj is the area of the face we have 


\Bi \ (c(xc,i''+^) - c{±c,i^)) 


r / - 


-f 

Pe J{k 




dc 

dx 


-A. 


Aw [vc] 
dc 


dx 


+ An [vc], 
dc 

d- An 

dy 


■ [vc]^^ + At [vc]^^ - ib [vc]^ J 


dr 


-A, 


dc 

dy 


dc 

dz 


- Ab 


dc 

dz 


dr. 


Denoting Cj(r) = c{'kj,r) and Vj(r) = v(xj,T) for j = n, s, e, w,t, b, c, by first order finite difference 
methods we have 


l^il (cc(P~^^) - Cc(P)) + / ^eVeCe - AwVwCw + ^nV„C„ - ^^VsC* + AtVtCt - AbVhCb dr 

£k + l 

f ^ f 'a ^6 ‘a ^A ^Ti 'a "a ‘a I t ^ ^ \ 

= / ( 22 ) 

Pe y ox ox oy oy oz oz J 

where 5x^ 6y and 6z are the width, length and height of the control volume Bi. By virtue of using 
voxelized geometry, we know that 5x = 5y = Sz, Aj = (Sx)^ for all j = n^s^e^wA^b^ and \Bi\ = {5x)^ . 
By the implicit Euler method, we, therefore, have 


Sx (( 


:/c+l 


At 






^/c+1 ^/c+1 ^/c+1 ^/c+1 




v/c+1-/c+l 


Pe 5x 


( c ^+1 + c ^+1 + cy 1 + c ^+1 


6c^i), (23) 


where c^ = c(xj, for j = n, 5, e, re, t, 6, c. In the case that one of the faces of the control volume lies on 
a boundary, the appropriate boundary conditions must be discretized; for the inlet, outlet and solid (or 
symmetry) boundaries this is straightforward due to the Dirichlet and zero Neumann boundary condition 
imposed there via (HU and ( |16gp . Therefore, we omit the details for the discretization of the boundary 
conditions on the external boundary dCl. The appropriate discretization of the reactive boundary con¬ 
ditions, prescribed on the fluid-solid interface through (jl6ep and the corresponding description of the 
reaction kinetics, here either dm, (HU) or dm, is slightly more involved and deserves a more detailed 
discussion. 
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In a fully implicit and coupled discreti z ation the resulting discrete equations are nonlinear and the 
Newton method needs to be used (|Kellevl. 119951) . In a broad class of practically interesting problems 
we have considered to date, we have not faced very strong coupling between the dissolved and adsorbed 
concentrations. Therefore, a fully implicit and coupled discretization was not required and we have found 
that an operator splitting approach, or just a Picard linearization, has worked well. In this approach, 
the dissolved particle concentration is computed at t = and then the value is used to compute the 

deposited mass at the time t = Runge-Kutta methods, or other methods for numerically solving 

stiff ODEs, ar e also straightforward to implement, and may be the subject of future studies if required 
( Kellevl . [ 1995 !) . In order to illustrate the method used in Pore-Chem, we describe the discretization for 
the Langmuir isotherm, which is achieved as follows. Firstly m is substituted into (|16ep . which is then 
split into a Robin boundary condition and an ordinary differential equation: 




1 - 


mo 


55 = Dal ,e fl - 

dt \ ^00,i / 


- Da^m, 

— DaVm. 


X e f> 0. 


(24a) 

(24b) 


If DaV-g- 1 - DaV = 0, then either c = DaV = 0 or DaV = DaV = 0. In both of these cases, 

^oo,i 

by m, fi = 0 and so no reactions occur at the spatiotemporal point under consideration, in which 
case, by ()24ap . we have Vc • n = 0 and a zero Neumann boundary condition, which is straightforward to 

implement. Otherwise, if Da^ -h DaJ^ ^ > 0, assuming that c(x, t) is constant over the time period 

in question, namely t G and equal to c(x) at each spatial point, (j24bp may be integrated to 

give 


m(x,P+^) 


Dakc(x) - -Bexp (- (Dayc(x)mJ^^ + Da^) 


DaV c(x)m^b + Da^ 


X e Tj, i > 0. 


Here H is a constant of integration which may be evaluated at t = to give 

B = (Dayc(x) (1 - 7h{x, - Daym(x, P)) exp ((Dayc(x)?h“b + Da^) P) 

Upon substitution into (l25l) . we have 


(25) 


(26) 


m(x, 

Dayc(x,P) - (DaUc(x,t'=) (l 


- Da\^im{x,P)) exp (- (Paj, 
Dai,iC(x, t'=)m“h + Da^_i 


+ Da^,i) (At)) 

5 

(27) 


for X G fi, t > 0, where we have approximated c(x) by This is done to prevent nonlinear terms 

in unknown variables from appearing in the discretized version of (jl6ep . Discretization of the other two 
isotherms is implemented in a similar manner. 

Consequently, we may approximate the Robin boundary condition, (jl6ep . on the reactive face J^ij G Pi 
for j = s^t^b and i = 0,l,...A^ — 1 using finite difference methods fully implicitly as follows: 


2n ^ I 


+ 2n 


,/c + l 




m 
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where n = ±1 is the direction of the outward pointing normal. Using (|28p . the appropriate terms are 
assembled, along with (|23P minus the relevant diffusive flux term, for the finite volume on which the 
reactive surface lies into a matrix and vector where and is a vector 

of the dimensionless solutions c(x, t^+^) at the discretized points of the computational domain. Once 
all the terms for all the finite volumes within the domain have been assembled into A^~^^ and the 





















equation is solved using a biconjugate gradient stabilized method to give the updated 

fluid concentration, at each discrete spatial point inCtf, while the updated adsorbed concentration, 
is given by (|26]). Time is then updated, the next timestep considered, and we proceed in the usual 
manner until the final time point is reached. 

Numerical simulation of the equations is performed in Pore-Chem and a schematic, outlining the 
numerical algorithm used, is given in Figure [H In the following section, the dimensionless equations are 
solved, but we present results in dimensional quantities. 


4. Results 


We present illustrative results, using the numerical method outlined in Section [3j on two separate 
computational domains, a real geometry and a virtually generated geometry, for two applications where 
surface reactions are important. The first set of simulations are performed on a portion of palatine sand¬ 
stone, obtained using micro con ipnterized tomograp hy (/i-CT) by Frieder Enzmann at the Johannes 


Gutenberg University of Mainz (|Becker et al.l . l2Qll[ ). Surface reactions are highly important in many 


fields of the Earth sc iences, includ i ng cal cite growth, oxidation-reduction reactions, formation of biofilms. 


to name only a few (|Steefel et al.l . |2QQ5[) . The second set of results is performed on a computational do¬ 


main virtually constructed to be representative of a commercially available microfiltration functionalized 
membrane. The use of functionalized membranes is a promising method for removing contaminants from 
water, an d involves treat ing the pore walls of the membrane so that they adsorb certain microorganisms 
or drugs (|Ulbrichtl . I2QQ6I) . Such membranes have pore sizes on the sub-micron scale and the resolution 
provided by ja-CT imaging techniques is not high enough to give representative images, motivating the 
use of a virtually generated geometry. The two computational domains under consideration are shown 
in Eigure[2l 

Eor the numerical simulations presented, a dead-end setup is used, with a schematic illustrating 
the domain and boundary conditions shown in Eigure[3l Layers of pure water voxels are added at the 
inlet and at the outlet, as shown in Eigure[2]and illustrated in Eigure[3j to allow free flow to develop. 
This results in a total computational domain size of 200 x 200 x 300 voxels for the sandstone geometry, 
and 100 x 100 x 120 voxels for the membrane geometry. Eor the results presented one type of reactive 
solid-fluid interface is considered, so that N = and consequently we now drop the i subscript notation 
relating to the type of reactive boundary. 


4 . 1 . Fluid flow 

As described in Section [3l we are first required to solve for the fluid flow. We assume that the fluid 
generates no slip as it passes over the membrane, and we set the slip length, /3, to be zero so that, by 
there is zero normal and tangential velocity at the fluid-solid interface, F. We set inflow velocities to be 
typical for the application under consideration, and use Un = 1 mm/s for the membrane geometry and 

= 1.5 X 10 mm/s for the rock geometry. The parameters chosen for the inlet velocities, and the fluid 
density and viscosity, yield small Reynolds numbers; Re = 4.2 x 10“^ in the case of the rock geometry 
and Re = 7.83 x 10“^ for the membrane geometry. Therefore, the fluid flow in both computational 
domains is in a Stokes regime. Solving (|14ap and (iMbll . along with the boundary conditions (|16ap - 
(|16dp . numerically until steady-state is achieved, yields the solutions as reproduced in Figure [H Due to 
our assumption that the maximal possible number of adsorbed particles is sufficiently small to ignore the 
effects of geometry modification, these remain constant through time. Examination of Eigure|4] reveals 
the dependence of the local velocity field on the morphology of the computational domains. 


4 . 2 . Contaminant transport 

We now turn our attention to solving the reactive flow. Illustrative parameters are used, and we 
employ the Henry isotherm, given by Equation m, for the rock geometry, and the Langmuir isotherm, 
given by Equation (HU), for the membrane geometry. We choose cin arbitrarily to be 2 x 10^ number/mm^ 
and we set the initial concentration of fluid contaminant to be equal to the inflow boundary condition, 
so that co(x) = 1 for X G Uj, while the quantity of adsorbed contaminant is initially assumed to 
be zero, mo(x) = 0. Eurthermore, for the membrane geometry, we set the dimensionless maximal 
surface concentration of adsorbed contaminant, ffioo, to be 10“^, which equates to a dimensional value 
of moo = 0.014 number/mm^. Due to the form of the equations, the choice of cin does not influence 
the dimensionless system of equations describing the transport and reaction of the contaminant given by 
(|14cp along with the boundary conditions (HSi - (HHl), except for through the parameter moo- Therefore, 
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Figure 1: Flow chart to illustrate how numerical computation of the transport with surface reactions is implemented 
within Pore-Chem, with main features indicated only. The light red ovals indicate input files, while the yellow ovals 
indicate output files. The green and the blue rectangular boxes indicate steps involved for the flow and efficiency solvers 
respectively, while the dark red boxes indicate steps involved in both. The white diamonds indicate decision making steps. 
In the case that the dimensional system of equations is being solved, the appropriate dimensionless quantities of interest 
are replaced by the dimensional equivalent. 
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(a) Rock geometry 



(b) Membrane geometry 


Figure 2: The two computational domains under consideration, plotted in 3D on the left-hand side of the 
figure, and in 2D through a representative cross section on the right-hand side of the figure. Figure [(a)] 
sh ows the palatine san dstone geometry obtained through /i-CT, with a voxel size of 1.4 x 10“^ m, as used 
in iBecker et al.l (|2Qll[) . Figure [(1^ shows a virtually generated geometry which aims to reproduce the 
morphology of a commercially available microfiltration membra ne, with a voxel leng th of 7.03 x 10“^ m 
(3 s.f.) This geometry was created using the soft ware GeoDict (IGeoDict. 2012-20141) . and a comparison 


to experimentally evaluated quantities is made in Di Nicolo et al. I (l2015[l . 
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Figure 3: Schematic to illustrate the computational domain. The solid microfiltration membrane is show 
in grey while the water is shown in white. Voxels are added to the top and botto m of the domain , at th e 
inlet and outlet, to enable free-fiow to develop. Modified, with permission, from iDi Nicolo et al.l (|2Q15f) . 



(c) Velocity magnitude, membrane (d) Pressure, membrane 

Figure 4: Velocity magnitude, |(a)| and |(c)[ and pressure, |(b)| and |(d)| fields, measured in mm/s and kPa 
respectively for both the computational domains considered. Due to our assumption that the adsorption 
of the solute does not alter the geometry, these remain constant throughout the experimental time frame. 
The black boxes shows the outline for the entire computational domain, where we exclude the additional 
voxels embedded at the inlet and the outlet for 3D plotting purposes. 
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Figure 5: The retained concentration of contaminant, m, for the rock geometry at 1200 second time 
intervals where Da^ = 0.1, Da^^ = 0.001 and Pe = 2.0xl0“^. The fluid concentration remains effectively 
constant throughout the domain and so this is not plotted. 



Figure 6: Numerical results at 5 x 10 ^ s intervals for the fluid concentration, c, and the adsorbed 
concentration, m, of contaminant in the membrane geometry. 


increasing or decreasing cin for fixed rhoo purely scales the dimensional concentration (both fluid and 
adsorbed) by a constant factor. 

Figure [5] illustrates the concentration of the contaminant on the solid surface over time for the rock 
geometry, where we use Da^ = 0.1, Da^^ = 0.001 and Pe = 2.0 x 10“^. Due to the slow rate of reaction 
and small Peclet number, the quantity of solute adsorbed over this time period is not large enough to 
significantly alter the fluid concentration over the time period under consideration. Furthermore, as 
the rate of reaction is much smaller than the mass transport rate, the adsorbed concentration remains 
spatially homogeneous. 

In contrast, with the parameters Pe = 10 and Daa = Da^^ = 10 for the membrane geometry, we 
see significant spatial heterogeneity in both the dissolved and adsorbed concentrations, as illustrated in 
Figure [6l As time progresses, the fast rates of reaction result in a depletion of the dissolved concentration 
at the pore wall and a dependence of both the dissolved and adsorbed concentrations on the local 
membrane morphology. 


5. Conclusions 

We have presented an algorithm for solving solute transport at the pore-scale within a resolved 
porous medium, with reversible surface adsorption at the pore wall. A pore-scale description of reactive 


13 
























transport, as opposed to a description at the Darcy scale, allows for a very accurate representation of 
the processes of interest. The system of equations comprise the NS equations and a CD equation, with 
Robin boundary conditions coupled to an ODE accounting for the surface reactions. Assuming that each 
particle is sufficiently small in size not to alter the flow of the fluid within the computational domain, 
and that its reaction at the wall does not significantly alter the pore-scale geometry, there is a one¬ 
way coupling between the NS and CD systems of equations. Although, for simplicity, we consider one 
species of solute and examine only surface reactions, extension to several different species of solute with 
both volumetric and surface reactions is straight forward and implemented within our software package 
Pore-Chem. The algorithm presented employs a FV method, and in this paper we have particularly 
focused on the discretization method used to solve the reactive boundary conditions for the adsorption 
and desorption at the interface. Illustrative numerical results, using our software package Pore-Chem, 
are presented on two separate geometries. The first of these is a 3D /i-CT image of a piece of Palatine 
Sandstone rock, while the second geometry is virtually generated within GeoDict ( GeoDict. 2012-20141) 
to be representative of a commercially available functionalized membrane. The results demonstrate the 
potential of such a numerical package, with the ability to solve reactive transport directly on images 
and on virtually generated geometries, in further progressing the understanding of the interplay between 
the transport and reaction rates at the pore-scale. In a future publication, currently in preparation, we 
will investigate the influence of the computational domain morphology on the reaction dynamics, and 
investigate the effect of different parameter regimes on the numerical results and quantities of interest. 

The advantages of using a pore-scale description are multiple. Firstly, it allows us to simulate reactive 
transport over a range of different parameter regimes, and in particular, outside the applicability region 
of the equivalent upscaled model. In highly disordered media the Peclet number can significantly vary 
locally, which poses problems for asymptotic upscaling methods, but not for a pore-scale description. 
Secondly, different kinetic models for the reactions can be used without the need to re-perform the 
upscaling procedure. In the future we plan to extend the algorithm in order to solve coupled mnlti scale 
problems using t he heterogeneous multiscale method, in a similar manner to iBattiato et al.l (|2Qllh and 
llliev et ^ ( 20141) . Such a development will enable problems at larger spatial scales to be considered, 
which could aid further research into the influence of pore-scale processes in a number of highly interesting 
and important research applications. 
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